The Medusa of Spatial Sorting: 
3D Kinetic Alpha Complexes and Implementation 
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Abstract 

Motivated by an application in cell biology, we consider spa- 
tial sorting processes defined by particles moving from an 
initial to a final configuration. We describe an algorithm 
for constructing a cell complex in space-time, called the 
medusa, that measures topological properties of the sorting 
process. The algorithm requires an extension of the kinetic 
data structures framework from Delaunay triangulations to 
fixed-radius alpha complexes. We report on several tech- 
niques to accelerate the computation. 

Keywords. Computational geometry, Delaunay triangulations, al- 
pha complexes, kinetic data structures, spatial sorting, exact geo- 
metric computation, implementation, software experiments. 



1 Introduction 

Consider a finite set of particles or points in M 3 , moving in 
time along continuous trajectories. Interpreting these points 
as the centers of moving objects, we are interested in the 
topological changes the configuration undergoes. Our in- 
terest in this problem originates in a sorting process that 
segregates cells during zebrafish development, as studied by 
Heisenberg and Krens J8). The sorting process operates on 
intermixed configurations of cells in which different types 
have different physical properties. One example is a mix 
of two cell types, in which the cells of the first type have a 
strong preference for neighboring cells of the same type and 
a strong dislike of exposed boundary, while the cells of the 
second type have milder preferences and dislikes. The typi- 
cal outcome is that the cells of the first type form a ball-like 
shape that is engulfed by a spherical shell consisting of cells 
of the second type. 



* This research is partially supported by NSF under grant DBI-0820624, 
by ESF under the Research Network Programme, and by the Russian Gov- 
ernment under mega project 11.G34.31.0053. 

1 1ST Austria (Institute of Science and Technology Austria), Klosterneu- 
burg, Austria. 

tlST Austria (Institute of Science and Technology Austria), Kloster- 
neuburg, Austria, Departments of Computer Science and of Mathematics, 
Duke University, Durham, North Carolina, and Geomagic, Research Trian- 
gle Park, North Carolina. 



In an effort to formalize the sorting process and to make 
it amenable to detailed and objective measurements, Heisen- 
berg, Krens, and the authors of this paper introduced the re- 
stricted Voronoi medusa as a mathematical representation. It 
is a geometric body in 4-dimensional space-time obtained by 
stacking up restricted Voronoi regions in M 3 J5). At any mo- 
ment in time, the Voronoi region of a particle is intersected 
with a ball, and the resulting bodies are glued together to 
form a 4-dimensional structure. Applying persistent homol- 
ogy to the time function on this structure yields fine-grained 
information about the sorting process that is difficult to ob- 
serve directly. This paper complements this foundational 
work with a description of the computational aspects of the 
medusa construction. 

Results. |5] proves that the restricted Voronoi medusa has 
the same homotopy type as the medusa obtained by stacking 
up simplices in the corresponding alpha complex. The latter 
alpha medusa is combinatorial in nature, which has compu- 
tational advantages. Continuing in this direction, this paper 
makes three contributions: 

1 . We describe an algorithm that maintains a fixed-radius 
alpha complex for points moving on piecewise alge- 
braic trajectories in R 3 . The algorithm supports inser- 
tions and deletions of points and allows for piecewise 
algebraic trajectories. 

2. Maintaining an alpha complex, we construct the alpha 
medusa whose geometric and topological properties re- 
flect the events during the sorting process. 

3. We convert the kinetic algorithm and the medusa con- 
struction into robust and efficient software. Basing the 
implementation on the Cgal package for kinetic data 
structures by Daniel Russel ifTTl . it achieves correctness 
through the exact geometric computation paradigm. 

Part of Contribution 1 is an extension of previous work 
from kinetic Delaunay triangulations to kinetic alpha com- 
plexes, which is of independent interest. The requirement 
of correctly comparing algebraic numbers, without tolerance 
for inaccuracy or approximation in Contribution 3 seriously 



slows down the software, even for piecewise-linear trajecto- 
ries. To counteract, we introduce techniques that speed-up 
the computations without sacrificing their correctness. We 
evaluate the effectivity of these techniques experimentally. 

Outline. Section [2] explains background from computa- 
tional geometry and topology. Section [3] describes the ki- 
netic algorithm for fixed-radius alpha complexes. Section |4] 
explains the algorithm that constructs the medusa of a set of 
moving cells. Section[5]describes techniques to speed up the 
computations. Section|6]concludes the paper. 

2 Background 

We review the fundamental geometric data structures that are 
required in this work. Voronoi tessellations and Delaunay tri- 
angulations are treated in most computational geometry text- 
books, including |2][3l, and alpha complexes are described 
in 13][6). For the most part, the discussion focuses on the 
3-dimensional case. Most definitions and properties extend 
to higher dimensions as well as to the plane. We will exploit 
the latter fact when we draw illustrations. 

Cell complexes. We recall that a k-simplex is the convex 
hull of an affinely independent set of k + 1 points in some 
Euclidean space. A face is a simplex defined by a subset 
of the k + 1 points. It is proper if the subset is different 
from the set. Reversing the direction, we call the fc-simplex 
a coface of its face. We define a simplicial complex as a finite 
collection of simplices that is closed under the face relation, 
with the additional property that any two simplices in the 
collection are either disjoint or their intersection is a face of 
both. The boundary of a fc-simplex is the collection of its 
(k — l)-faces. The simplices of dimension 0, 1, 2, and 3 
are referred to as vertices, edges, triangles, and tetrahedra. 
The star of a fc-simplex is the set of simplices that contain 
the fc-simplex as a face. Noting that the star is in general not 
closed under the face relation, we define the closed star as 
the set of all simplices in the star and of all faces of these 
simplices. It is the smallest simplicial complex that contains 
the star. Finally, if a is a fc-simplex and u is a point that does 
not lie in the fc-plane of the simplex, then the join, denoted 
as u * a, is the (fc + l)-simplex that is the convex hull of u 
and the vertices of a. 

It is convenient to also introduce an abstract counterpart to 
the above geometric concept of a simplicial complex. Specif- 
ically, an abstract simplicial complex consists of a finite 
set of (abstract) elements and a collection of subsets that is 
closed under the subset relation. We may map each element 
to a point in some Euclidean space, and each subset to the 
convex hull of the points that correspond to its elements. If 
the dimension of the space is sufficiently high and the points 
are well chosen, this is a simplicial complex, which we re- 
fer to as a geometric realization of the abstract simplicial 
complex. Here is an example of this construction. Consider 



a finite set, X, of possibly overlapping bodies, and define 
the nerve as the collection of subsets of X with non-empty 
common intersection. We note that the nerve is an abstract 
simplicial complex. Indeed, the bodies are the elements, and 
if A C X is a set in the nerve, then every subset of A is 
also in the nerve. A useful result is the Nerve Theorem J4], 
which states that if the bodies in X are convex then every 
geometric realization of the nerve has the same homotopy 
type as the union of the bodies. Intuitively, this means that 
one can be transformed into the other by continuous transfor- 
mations like bending, shrinking, and expanding, but without 
gluing and cutting. Most complexes in this paper will be 
constructed using the nerve operation. 

We will have occasion to also use complexes that are more 
general than simplicial complexes. Instead of fc-simplices, 
they contain k-cells, which are homeomorphic to the fc- 
dimensional unit ball. The boundary of a fc-cell is the home- 
omorphic image of the (fc — l)-sphere that bounds the fc-ball. 
We define a cell complex as a finite collection of cells with 
pairwise disjoint interiors such that the boundary of each fc- 
cell is a union of (fc — l)-cells in the complex. 

Voronoi tessellations and Delaunay complexes. Con- 
sider now a finite set of points, U, in M 3 . The Voronoi region 
of a point u in U is the set of points x G R 3 that have u as 
the closest point in U : 

) = {x G M 3 | \\x - u\\ < \\x - v\\,Vv G U}. 

Note that votjj(u) is convex. We usually drop U from the 
notation. The Voronoi tessellation of U is the set of Voronoi 
regions of its points. While it is not a cell complex formally, 
we get a cell complex if we add the common intersections of 
Voronoi regions as lower-dimensional cells to the set. If the 
points in U are in general position, by which we mean that 
no four lie in a common plane and no five lie on a common 
sphere, then the Voronoi regions intersect in a rather pre- 
dictable pattern. Specifically, the intersection of any two is 
either empty or a (2-dimensional) polygon, the intersection 
of any three is either empty or a (1 -dimensional) edge, and 
the intersection any four is either empty or a (0-dimensional) 
point. Furthermore, the intersection of five or more Voronoi 
regions is necessarily empty. 

We get the dual Delaunay complex if we replace each non- 
empty intersection of Voronoi regions by the convex hull of 
the points that generate the Voronoi regions containing the 
intersection. Equivalently, we may define the Delaunay com- 
plex as the set of convex hulls of subsets of points that have 
the empty sphere property. Specifically, this means that there 
exists a sphere that passes through the points of the subset 
and all other points in U lie strictly outside this sphere. We 
note that the center of this sphere belongs to the intersection 
of the corresponding Voronoi regions. Assuming general po- 
sition, the Delaunay complex is a simplicial complex, which 
is generally referred to as the Delaunay triangulation. It is a 
geometric realization of the nerve of the Voronoi tessellation. 
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Restricted Voronoi tessellations and alpha complexes. 

Fixing a positive radius, ao, we define the restriction of a 
Voronoi region to be its intersection with the closed ball of 
radius ao centered at the generating point: 

res[/(u) = {x G votu(u) | ||x — u\\ < ao}. 

Again, we usually drop U from the notation. The restricted 
Voronoi tessellation of U is the set of restricted Voronoi re- 
gions of its points. In contrast to the unrestricted case, each 
restricted Voronoi region is bounded, and the tessellation 
covers only the union of balls and not the entire space. 

As before, we assume general position so we can dualize 
by geometrically realizing the nerve. The resulting simpli- 
cial complex is called the alpha complex. Since res (it) C 
vor(u), for each point u in U, the alpha complex is a sub- 
complex of the Delaunay triangulation. Next, we derive an 
equivalent condition for a Delaunay simplex to lie in the al- 
pha complex which is more suitable for computations. Each 
fc-simplex in the Delaunay triangulation has a unique cir- 
cumscribed (fc — l)-sphere in its supporting fc-plane. We 
call its center the circumcenter, its radius the circumradius, 
the ball in R 3 with this center and this radius the circumball, 
and the sphere that bounds the circumball the circumsphere 
of the fc-simplex. Note that the circumsphere is the smallest 
sphere that passes through the vertices of the fc-simplex. We 
call the fc-simplex short if its circumradius is smaller than 
or equal to ao- Finally, we call the fc-simplex Gabriel if its 
circumball has no point of U in its interior. 

Short&Gabriel Lemma. A simplex in the Delaunay 
triangulation of U belongs to the alpha complex, for radius 
ao, iff it is short and Gabriel, or it is the face of another 
Delaunay simplex that is short and Gabriel. 

The face of a short simplex is necessarily short, but the face 
of a Gabriel simplex is not necessarily Gabriel. It follows 
that all simplices in the alpha complex are short, but not all 
simplices need to be Gabriel. Also note that a tetrahedron 
is in the Delaunay triangulation iff it is Gabriel; therefore, 
it is in the alpha complex iff it is short. In our application, 
we use the restricted Voronoi tessellation to model a set of 
biological cells for which the positions of their nuclei are 
known. Indeed, a cell tends to minimize its surface area and 
usually does not grow larger than a certain size. Therefore, a 
restricted Voronoi region appears to be a good approximation 
of the actual cell shape and is still simple enough for our 
computational purposes. 

3 Kinetic Alpha Complexes 

In this section, we describe the algorithm that maintains the 
alpha complex for a fixed radius ao > 0. We pay particular 
attention to the certificates that govern the sequence of oper- 
ations needed to preserve the correctness of the structure at 
all times. 



The kinetic framework. The input to our algorithm is a 
finite set of trajectories, each a continuous map u : [a, b] —> 
R 3 with < a < b < 1. For simplicity, we assume it to be 
piecewise linear, with a = t n < t\ < . . . < tf. = b such that 
there are points a,j , bj £ R 3 for which u(t) = (1 — t)aj + tbj 
for tj < t < tj+i. In other words, we can write u(t) = 
(/i(*)> fz(t)i fs(t)) sucn that between tj and tj+i, each fi is 
a polynomial of degree 1. We call to, ti, • • • , tk the bending 
events of the trajectory. Furthermore, we assume that the 
trajectories do not meet each other, that is, u(t) ^ u'(f) for 
all u, u' G U and all t for which both trajectories are defined. 

Our task is to maintain a data structure that goes from an 
initial configuration, at time t = 0, to the final configuration, 
at time t = 1. For that, the data structure is constructed at 
time t = 0, and maintained through a sequence of update op- 
erations until the final configuration is reached. It is assumed 
that the number of updates is finite, and we call the time of 
an update an event. Events are detected by defining suitable 
certificate functions, also referred to as certificates. At any 
moment t different from any event, we have a collection of 
active certificates, all being non-zero at t. Importantly, they 
guarantee that as long as no certificate changes its sign, our 
data structure remains structurally unchanged. To detect the 
next event, the algorithm then finds the smallest root of any 
active certificate that is greater than t. It handles the event 
by updating the data structure and the collection of active 
certificates. Throughout this paper, we make the simplifying 
assumption that all events are distinct, that is, no two events 
happen at the same moment in time; see also Section|6] 

Maintaining the Delaunay triangulation. Since we need 
it later, we begin by reviewing the kinetic algorithm for 3- 
dimensional Delaunay triangulations described in lfIT'1. Be- 
sides changes brought about by insertions and deletions of 
points, and switches to new trajectory segments, there are 
only two configurations that trigger a structural change in 
the triangulation: 

• five points of U lie on a common sphere, and the open 
ball bounded by this sphere contains no points of U\ 

• four points of U lie on a common plane, and one of 
the open half-spaces bounded by this plane contains no 
points of U. 

We call each such configuration a degeneracy. Consistent 
with the above assumption of distinct events, we assume that 
at every moment of time there is only one degeneracy, and 
that each degeneracy lasts only for a single moment. In other 
words, we can find a small open interval in time during which 
the given degeneracy exists at a single point in time, and it 
is the only degeneracy that occurs during this interval. We 
can therefore study the effect of the degeneracy by consider- 
ing the non-degenerate local configurations right before and 
right after the degeneracy. Consider for example a degener- 
acy of the first type, which involves five points. Right before 
the degeneracy, the five points span two Delaunay tetrahedra 



3 



with a common triangle, and right after the degeneracy they 
span three tetrahedra so that each pair shares a triangle and 
all three share an edge. Of course, it can also be the other 
way round. Importantly, we can transform one configuration 
to the other by flipping. In this particular case, we substitute 
three for two or two for three tetrahedra, calling the opera- 
tion a 2-3-flip; see Figure [T] To avoid a case analysis, we 




Figure 1: Illustration of a 2-3-flip that alters the triangulation of 
a triangular double pyramid. On the left, the five points span two 
tetrahedra meeting in a triangle. After the flip, the triangle is re- 
placed by an edge and the three incident triangles that connect the 
edge to the remaining three points. 

represent the triangulation using a vertex at infinity that is 
joined to every simplex in the boundary of the convex hull 
of U. Effectively, we embed the triangulation on a 3-sphere. 
This way, we can add the point at infinity to the set of four 
points forming a degeneracy of the second type, thus getting 
a degeneracy of the first type, which is handled by a 2-3-flip, 
as described above. 

Flip events. The transition of the Delaunay triangulation 
across degenerate configurations is controlled by two cer- 
tificate functions. Let u 1 ,u a ,'U a , , u 4 ,'U. 5 be the five trajec- 
tories of the points that span two tetrahedra sharing a trian- 
gle or three tetrahedra sharing a common edge, as in Figure 
Q] If one of the trajectories belongs to the infinite vertex 
then we reorder them such that this trajectory is u 5 . Let 
= /K*); be the coordinate functions of 

the finite points, and recall that the squared norm of this point 
is the sum of the squares of its three coordinates. If all five 
points are finite, we create the certificate 
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which is a univariate polynomial in t that is zero iff the five 
points are co-spherical. Assuming the coordinate functions 
are linear, the degree of the polynomial is 5. If the fifth point 
is at infinity, we create the certificate 

" i flit) fiit) fiit) - 

det i fm fiit) fi(t) 

i fiit) fiit) fS(t) ' (2) 

i fiit) fiit) fi(t) 



which is zero iff the four point are coplanar. We call the poly- 
nomials in (1112b flip certificates and their roots flip events. 

After having constructed the initial certificates, at time 
t = 0, the algorithm finds the first positive flip event. It 
then performs a 2-3-flip, creating certificates for the (new) 
simplices inside the double pyramid, and updating the certifi- 
cates of the simplices in the boundary of the double pyramid. 
The updating is necessary because the star of each boundary 
simplex changes during the flip. After these steps, both data 
structure and certificates are again valid, and the iteration 
continues with the next flip event. 

Radius events. Next, we extend the kinetic algorithm from 
Delaunay triangulations to alpha complexes. As before, we 
use a fixed radius ao > 0. We represent the alpha complex 
by equipping each Delaunay simplex with a flag that indi- 
cates whether or not the simplex belongs to the alpha com- 
plex. To construct these flags at time t — 0, we check every 
Delaunay simplex for being short and for being Gabriel. Fol- 
lowing the Short&Gabriel Lemma in Section [2] we add all 
Delaunay simplices that are short and Gabriel, as well as all 
their faces, to the alpha complex. To maintain the flags, we 
construct a certificate for each edge, triangle, and tetrahedron 
whose roots are the times when the circumradius of the sim- 
plex equals ao- To simplify the discussion, we assume the 
generic case in which the circumradius changes from strictly 
smaller to strictly larger than ao, or vice versa. We call these 
functions radius certificates and their roots radius events. 
Whenever a Delaunay simplex is inserted or deleted, the al- 
gorithm also creates or removes the corresponding radius 
certificate. The certificate of an edge compares the length 
to 2«o> an d taking squares, we get a polynomial of degree 
2. The radius certificates of a triangle and a tetrahedron are 
more complicated, but can be derived from suitable minors 
of the matrix that defines the circumsphere of the simplex; 
see lfl3l for the formula in the tetrahedral case. We will dis- 
cuss the triangle case in Section 

After initializing the alpha complex and the certificates, 
the algorithm looks for the next event. If this is a flip event, 
we proceed as described above. In addition, we update the 
flags that identify the alpha complex as a subcomplex of the 
Delaunay triangulation. Because all tetrahedra involved in 
the 2-3-flip have the same circumsphere, they are either all 
short or all non-short. If they are short, all new Delaunay 
simplices are added to the alpha complex, and otherwise, 
none of them is added to the alpha complex. Second, con- 
sider the case in which the next event is a radius event. Let a 
be the conesponding Delaunay simplex. If a goes from non- 
short to short, then its proper cofaces are necessarily non- 
short. We check whether a is also Gabriel, noting that this 
is always the case if a is a tetrahedron. If so, a is added to 
the alpha complex together will all its faces. On the other 
hand, if a goes from short to non-short, then its proper faces 
are necessarily short. We remove a from the alpha complex, 
unless is was not in the complex even before the event. If the 
event causes the deletion of a from the alpha complex, then 
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this may have consequences for its faces. In particular, if a 
was the last proper coface of r in the alpha complex, and r 
is not Gabriel, then r is also deleted from the alpha complex. 
Afterwards, the algorithm continues with the next event. 

Redundancy of Gabriel events. Perhaps surprisingly, flip 
and radius events suffice to maintain the alpha complex. Flip 
certificates monitor when tetrahedra become non-Delaunay, 
and radius events monitor when simplices become short or 
non-short. We do not need certificates that monitor when 
simplices become Gabriel or non-Gabriel. To understand 
why such certificates are not necessary, we call a time t 
G-critical for a simplex a, if at that time, a changes from 
Gabriel to non-Gabriel, or vice versa. 

G-CRITICALITY Lemma. Let t be a G-critical time for a 
short Delaunay edge or triangle. Then this edge or triangle 
has a proper coface that is in the alpha complex at time t. 

PROOF. Denote the edge or triangle by a and consider its 
circumball, B a , at time t. No point of U lies in the interior 
of B a , but there is a point u on the bounding sphere that is 
not a vertex of a. The join u * a is another simplex in the 
Delaunay triangulation, and it is a proper coface of a. It has 
the same circumball as a, which implies that u * a is short 
and Gabriel and therefore belongs to the alpha complex at 
time t. ED 

The lemma implies that when a short edge or triangle 
changes its Gabriel status, it is a face of a simplex in the 
alpha complex. The status change has therefore no impact 
on its membership in the alpha complex. 

Other events. For later reference, we briefly mention the 
remaining types of events supported by our algorithm. First, 
we consider a bending event, at which a trajectory starts 
a new segment. Such an event leaves the alpha complex 
unchanged, but all flip and radius certificates that involve 
the coordinates of the corresponding vertex are recomputed. 
These are the certificates associated to the simplices in the 
closed star of the vertex. 

Second, we allow for insertions and deletions of points. 
The two operations are mostly symmetric, and we only dis- 
cuss the insertion of a point u. We add u to the Delaunay 
triangulation by identifying all tetrahedra whose circumballs 
contain u, referring to their union as the conflict region of 
u. Since these tetrahedra no longer satisfy the empty sphere 
criterion, we remove them from the Delaunay triangulation, 
together with their faces in the interior of the conflict region. 
Next, we add u and connect it to all simplices in the bound- 
ary of the conflict region. After the operation, these sim- 
plices form the boundary of the closed star of it. Finally, the 
simplices in the closed star are checked for being in the alpha 
complex, and their certificates are created or updated. 



4 Medusa Construction 

In this section, we focus on the construction of space-time 
complexes that give a static 4-dimensional representation of 
a dynamic spatial process. After a brief review of the medusa 
concept, we discuss its construction using kinetic algorithms. 

Restricted Voronoi and alpha medusa. Let U be a finite 
set of trajectories, and write U(i) C R 3 for the set of points 
obtained by evaluating all trajectories at time t. Of course, 
we evaluate only those trajectories that are defined at t. The 
restricted Voronoi tessellation of U(£) is constructed as de- 
scribed in Section|2] Letting u : [a, b] — > R 3 be a particular 
trajectory, we get the region of u(i) within the tessellation 
of U(t) for each time a < t < b. Piling up these regions in 
R 4 and taking the closure, we call the result the stack of the 
trajectory u: 

stack (u) = clos I I^J [resu(t)(u(i)) x t] J . 

\a<t<b J 

Taking the closure is not redundant because the region of 
a trajectory can change discontinuously when inserting an- 
other trajectories, causing a locally open pile. The set of 
stacks of trajectories in U is the restricted Voronoi medusa, 
which we denote as 1Z = 1Z(U). While it is the preferred 
representation of a spatial sorting process, we introduce a 
dual structure that is easier to compute and gives the same 
measurements. As a first approximation of this construction, 
we exploit the duality between restricted Voronoi diagrams 
and alpha complexes and pile up the simplices of the latter 
to form cells in R 4 ; see Figure|2] This construction has two 




\ 

Figure 2: The triangle spanned by points on the three trajectories 
sweeps out a triangular prism while it belongs to the alpha complex 
between the lower and the upper planes. 

shortcomings. First, the cells do not define a cell complex, 
but this can be corrected by adding (geometrically degener- 
ate) cells that represent flip, insertion, and deletion events. 
Second, the walls are redundant and are better contracted. 
The contraction produces simplices, but they are connected 
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to each other in more general ways than allowed in a simpli- 
cial complex. 

The main difficulty in defining the alpha medusa formally 
is a lack of regularity of the stacks in 1Z. While every re- 
stricted Voronoi region is convex, the stack of such regions 
is generally not convex, and the intersection of two stacks is 
generally not connected. Taking the nerve of the set of stacks 
would thus produce an abstract complex that does not respect 
the Nerve Theorem. To finesse these difficulties, we take the 
multi-nerve as follows. If k + 1 stacks intersect in I different 
connected components, then the multi-nerve has I copies of 
the abstract fc-simplex, each representing one component of 
the intersection. Moreover, each abstract simplex is associ- 
ated with the lifetime interval, [ti,^], of the corresponding 
connected component. Because there is at most one com- 
ponent at every moment in time, we have disjoint intervals 
for different copies of the same simplex. The case t\ = ti 
is allowed and represents connectivity at a single moment 
in time. The result of this multi-nerve construction is what 
we call the alpha medusa and denote as A = A(XS). In J5J 
Lemma B], it is shown that A is homotopy equivalent to 1Z 
as well as the space covered by the prismatic cells swept out 
by simplices of the alpha complex. 

Construction. A simplex with lifetime interval [t\ , t?\ is 
active for all t\ < t < t?, dead finished for all t% < t. We con- 
struct the alpha medusa by running the kinetic alpha complex 
algorithm and updating the medusa at each event. We main- 
tain two lists of simplices, called the output list and the active 
list, preserving the following property at all times: 

Invariant. After the algorithm has handled an event at 
time t, the output list consists of all simplices finished at t, 
and the active list consists of all simplices active at t. 

At the beginning, all simplices in the alpha complex are put 
into the active list. Whenever a simplex is added to the alpha 
complex at a time t, we put it into the active list, and when 
it is removed from the alpha complex at a later time t', we 
remove it from the active list and put it into the output list 
with lifetime interval [t, t']. 

There are additional steps to be taken during special types 
of events. We discuss flip events now and the more compli- 
cated insertion and deletion events later. Consider a 2-3-flip 
at time t. We have two tetrahedra before and three after the 
flip, or vice versa. At time t, they all have the same circum- 
ball, so either all or none of them belong to the alpha com- 
plex. In the former case, the five corresponding restricted 
Voronoi regions meet in the center of the common circum- 
ball. It follows that the dual 4-simplex belongs to the multi- 
nerve, with lifetime interval [t, t\. We thus add the 4-simplex, 
connecting it its five boundary tetrahedra, which are the ones 
involved in the flip. Indeed, the 4-simplex fills the void be- 
tween the five tetrahedra, which could not be filled by gluing 
the tetrahedra to each other as they are not face-to-face. 



Insertions and deletions. Finally, we consider the case in 
which a trajectory, u', is inserted at time t. The case of a 
deletion is symmetric and omitted. Let u' = u'(t) be the ver- 
tex inserted into the alpha complex. As described in Section 
[3] the insertion of u' into the Delaunay triangulation is syn- 
onymous with the substitution of the star for the simplices in 
the conflict region of v! . Any subset of these simplices may 
belong to the alpha complex. In addition to finding this sub- 
set, we need to construct simplices that fill the gap between 
the deleted and the inserted alpha complex simplices in the 
medusa. This is similar to a flip, in which the void formed by 
the five tetrahedra is filled by a 4-simplex. To describe a gen- 
eral solution to this problem, let U be a finite set of points, 
v! a point not in U, and write U' = U U {u'} for the set 
including the new point. We call each rest; (it) an old region, 
each res;/' (it) a new region, and resu> (u') the region of u'. 

Insertion Lemma. Let I be a non-empty common in- 
tersection of old regions. Then, / intersects the region of 
v! iff the common intersection of the corresponding new re- 
gions is either empty or it also intersects the region of vl . 

Proof. For the forward direction, assume that / intersects 
the region of v! , and let x be an arbitrary point in that in- 
tersection. Note that \\x — u'\\ < d(x), where d(x) is the 
distance to the closest points in U, In the first case, we have 
\\y — u'\\ < d(y), for all y £ I, which implies that the cor- 
responding new regions do not intersect. In the second case, 
we have \\y — u'\\ > d(y), for some y, Then there exists a 
point y' G I such that \\y' — u'\\ = d(y'), which implies 
that y' belongs to the intersection of the new regions with 
the region of v! . 

For the backward direction, if the new regions intersect the 
region of u', so do the old regions. In the remaining case, the 
new regions do not intersect. Let a; be a point in the intersec- 
tion of the old regions. Since x is not in the intersection of 
the new regions, it must be in the region of v! . EQ 

We use the Insertion Lemma to derive an algorithm that 
updates the alpha medusa upon the insertion of a point v! at 
time t. Similar to the case of the Delaunay triangulation, we 
identify the conflict region and remove all simplices in its in- 
terior. In addition, if a removed fc-simplex a is in the alpha 
complex just before t, we add the join, u' * a, to the alpha 
medusa with lifetime interval [t,t]. We do this because the 
k + 1 old regions represented by a have a non-empty com- 
mon intersection, but the corresponding new regions do not. 
In the second step, we insert the star of u' into the Delaunay 
triangulation, check for every inserted simplex whether it be- 
longs to the alpha complex, and if it does then we add it to 
the active list. Clearly, every inserted simplex is of the form 
u' * a, with a on the boundary of the conflict region. 

5 Implementation and Experiments 

In this section, we turn to implementation issues. In par- 
ticular, we discuss how to implement the algorithm in a ro- 
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bust manner, we study the effect of practical choices, and we 
present experimental results obtained with our software. 

Robust computation. Recall the basic structure of a ki- 
netic data structure as explained in Section [3] it consists of 
certificate functions, which are polynomials in t, and each 
step advances the state to the smallest root larger than the 
current time. To maintain the certificates and to advance to 
the next event, the algorithm computes and compares real 
roots of univariate polynomials. These roots are algebraic 
numbers — irrational in general — which makes the com- 
putations non-trivial. The naive solution of approximating 
these roots by inexact floating-point numbers can have un- 
predictable effects. It is not true that the outcome is just 
slightly wrong, e.g. by switching the order of events that hap- 
pen almost simultaneously, but the incorrect order can lead to 
inconsistent configurations, causing program crashes, non- 
termination, and non-sensical results. This problem is well- 
known in geometric contexts ifTOl and several approaches 
have been proposed. We follow the exact geometric com- 
putation (EGC) paradigm, popularized by Chee Yap |fl5l . It 
suggests that the basic primitives be mathematically correct, 
so that an algorithm using these primitives is in the position 
to compute provably correct results. Translated to our situa- 
tion, we require that the events of our process are handled in 
the mathematically correct order. The price we pay for this 
interpretation of robustness is the burden to compute with 
algebraic numbers. 

We implement our algorithm using the CGAL librarjQ, 
which is designed in the spirit of the EGC paradigm. An- 
other aspect of Cgal is its generic programming approach: 
algorithms access underlying data structures and primitives 
through a well-defined interface, so that these layers can 
be easily replaced with alternative implementations. More 
specifically, we make use of the kinetic data structures pack- 
age lfl2l . which provides an EGC implementation of kinetic 
Delaunay triangulations in two and three dimensions. In- 
ternally, the package contains an algebraic kernel, providing 
the low-level functionality needed to handle roots of polyno- 
mials, and a combinatorial layer, maintaining the data struc- 
ture and the certificates over time. As mentioned earlier, we 
have extended the combinatorial layer to maintaining an al- 
pha complex. 

Experimental set-up. We use datasets obtained with the 
CompuCell3D software^, which allows for the simulation 
of a 3-dimensional cell segregation process using a Monte- 
Carlo algorithm for energy minimization; see the companion 
paper [5 | for more details. We focus on simulated as opposed 
to observed data because they offer a better control of the in- 
put size and the direct accessibility of the cell trajectories. 
In our particular example, the cells are colored blue or red, 
each color with probability one half, and the parameters of 

'Computational Geometry Algorithms Library, www . cgal . org. 
2 www . compucell3d . org/. 



the simulation are chosen so that the blue cells eventually en- 
gulf the red ones; see Figure[3]for an illustration. We created 




Figure 3: The restricted Voronoi tessellation at four moments in 
time. At the beginning, the cells form a cubical grid (upper-left). 
The cells move toward the center of the available space, and the 
blue cells begin to engulf the red cells (upper-right and lower-left), 
allowing for satellites while this happens. Finally, the blue cells 
form a sphere surrounding a ball of red cells (lower- right). 

datasets for several input sizes. In all cases, each trajectory 
represents the path of a cell nucleus which exists through- 
out the entire process. Hence, no new cells are ever inserted 
after the start of the process, and no old cells are deleted be- 
fore the end of the process. The trajectories follow a global 
rhythm in which each trajectory starts a new segment at each 
value in a common sequence of bending events. In between 
two bending events, each trajectory is linear. All experiments 
are performed on a Intel Core 2 Dual CPU clocked with 2.4 
GHz each, with 3 MB of cache size, and 4 GB of total mem- 
ory. The code runs under Debian Squeeze, compiled with 
gcc-4.4.5 and Cgal version 3.9. 



Input #s 




Time in sec 


#EVENTS 


traj bends 


Del 


Alpha 


Medusa 


flips 


rad 


20 20 


7 


740 


743 


512 


631 


20 40 


29 


1,550 


1,553 


1,011 


1,335 


20 80 


81 


3,205 


3,203 


2,019 


2,503 


20 160 


188 


6,473 


6,484 


3,978 


4,506 


10 40 


7 


487 


491 


369 


554 


20 40 


29 


1,549 


1,556 


1,011 


1,335 


40 40 


79 


3,975 


3,985 


2,874 


2,171 


80 40 


229 


9,897 


9,904 


7,856 


4,977 


160 40 


495 


21,516 


21,741 


17,667 


5,998 



Table 1 : Columns from left to right: the number of trajectories and 
bending events per trajectory, the time to maintain the Delaunay 
triangulation and the alpha complex, the time to compute the alpha 
medusa (which includes the maintenance of the alpha complex), 
and the number of flip and radius events. 

In a first test, we compare the running times for maintain- 



7 



ing the Delaunay triangulation and the alpha complex, and 
for constructing the alpha medusa, all for the same datasets; 
see Table Q] We observe that the overhead of computing the 
medusa is negligible. For this reason, we concentrate on the 
maintenance of the alpha complex. Comparing the third and 
fourth columns of the table, we see that the radius events 
slow down the algorithm by more than a magnitude, in spite 
of the fact that their number is not much larger than the num- 
ber of flip events. In the remainder of this section, we explain 
improvements of our algorithm aimed at reducing the perfor- 
mance gap between Delaunay and alpha complexes. 

Number of certificates. The bottleneck is the construction 
of radius certificates and the computation of their real roots. 
Recall that in our original formulation, we maintain a radius 
certificate for each edge, triangle, and tetrahedron. Our first 
optimization is based on the observation that many of these 
certificates are not necessary: if a simplex is short, then all 
its faces are short, and if a simplex is non-short, then all its 
cofaces are non-short. 

Optimization 1. Whenever a triangle or tetrahedron 
becomes short, we remove the radius certificates of its proper 
faces, and when an edge or triangle becomes non-short, we 
remove the radius certificates of its proper cofaces. 

Of course, this implies that we sometimes have to construct 
certificates that would otherwise still exist. For example, we 
construct the certificate of a triangle at the time its third edge 
becomes short. On the other hand, we avoid unnecessary cer- 
tificates, for instance the certificates of the boundary edges 
of a triangle that stays short for the whole simulation. As 
we see in Table [2] the strategy saves time in practice. We 
observe that the constructions of radius certificates and the 
running time both decrease roughly by a factor of two. 



Input #s 


Time in sec 


#Certificates 


traj 


bends 


before 


after 


before 


after 


20 


20 


740 


361 


20,211 


10,262 


20 


40 


1,550 


770 


40,897 


20,622 


20 


80 


3,205 


1,579 


82,287 


41,037 


20 


160 


6,473 


3,142 


163,511 


80,489 


10 


40 


487 


248 


12,932 


6,324 


20 


40 


1,549 


770 


40,897 


20,622 


40 


40 


3,975 


1,892 


105,754 


54,426 


80 


40 


9,897 


4,882 


259,848 


139,816 


160 


40 


21,516 


10,181 


566,589 


303,065 



Table 2: Third and fourth columns: the time to maintain the alpha 
complex before and after Optimization 1. Fifth and sixth columns: 
the number of radius certificates before and after Optimization 1 . 



Degree. We turn to the computation of certificates. As- 
suming piecewise-linear trajectories, the radius certificate of 
an edge is a polynomial of degree 2; compare with Section 
|3] There is a standard construction of a radius certificate of 
a tetrahedron, which is a polynomial of degree 8; see iTPJl . 
Our interest lies in the remaining triangle case. The current 



Cgal implementation computes the squared circumradius 
of a triangle A in R 3 with an expression of the form 

2 nmiij + num^ + nurn^ 



where den is the determinant of the matrix in (0, and num x , 
nuiriy, num 2 are expressions formed by minors of this ma- 
trix. The corresponding certificate, 

2 2 2 2 2 

mim x + miuiy + num z — 4a dcn , 

is a polynomial whose degree is 10, which is higher than the 
degree for the tetrahedron. We replace ([3]) by a simpler ex- 
pression. Writing u, v, w for the three vertices of the triangle, 
the circumradius can also be written as 

\\u-v\\ ■ \\u-w\\ ■ \\v-w\\ 

7a = — ^TiT ^ 7 Til ' W 

2||(u — w) x [v — w)\\ 

a formula that is straightforward to derive using elementary 
matrix calculus; see also Wikipedia fl4| . 

Optimization 2. Monitor the radius of a triangle using 
the following certificate function: 

\\u — w\\ 2 \\u — w\\ 2 \\v — w\\ 2 — 4ao||(u — w) x (v — w)\\ 2 . 

The degree of this certificate is 6. We see the effect of this 
improvement in Table|3] The running time improves by more 
than a factor of two. 



Input #s 


Time in sec 


traj 


bends 


deg 10 


deg 6 


20 


20 


361 


151 


20 


40 


770 


342 


20 


80 


1,579 


734 


20 


160 


3,142 


1,515 


10 


40 


248 


105 


20 


40 


770 


344 


40 


40 


1,892 


912 


80 


40 


4,882 


2,374 


160 


40 


10,181 


5,157 



Table 3: Timings for maintaining the alpha complex using a degree 
10 versus a degree 6 certificate function for monitoring the circum- 
radii of triangles. 



Algebraic kernel. As already mentioned, the Cgal pack- 
age for kinetic data structures contains an internal algebraic 
kernel, which, among other things, is used to isolate the 
roots of polynomials and sort them in the event queue. By 
the generic design of the package, the combinatorial layer 
communicates with the kernel via a small and well-defined 
interface, which makes it possible to replace the algebraic 
kernel with a different implementation. In recent years, a 
mature and generic algebraic kernel for geometric compu- 
tations has been developed [fl]. It has been integrated into 
CGAL and is available since version 3.7 under the name 
Algebraic_kernel_d. We refer to it as the ak_d kernel. 
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Input #s 




Time in sec 




traj 


bends 


kds 




filter 


cache 


20 


20 


i ^ i 

i j i 


84 


72 


4 / 


20 


40 


342 


176 


152 


98 


20 


80 


734 


364 


310 


198 


20 


160 


1,515 


731 


622 


390 


10 


40 


105 


55 


49 


30 


20 


40 


344 


177 


152 


100 


40 


40 


912 


458 


392 


256 


80 


40 


2,374 


1,180 


1,024 


689 


160 


40 


5,157 


2,566 


2,256 


1,481 



Table 4: Timings for maintaining an alpha complex using the kds 
kernel, the ak_d kernel, the latter with Descartes filtering, and in 
addition with enabled cache. 

Both the internal kds and the ak_d kernels use subdivision 
methods for root isolation, but they differ in the strategy for 
detecting empty intervals and isolating intervals. The kds 
kernel uses Sturm theory |[T6l §7], while the ak_d kernel is 
based on Descartes' rule of sign Q, which leads to a better 
performance in practice; see (§1 for a comparison of vari- 
ous root solvers. The difference between the third and fourth 
columns in Table [4] shows that exchanging the kernel yields 
another improvement of roughly a factor of two. 

Filter and cache. We get further optimizations by exploit- 
ing the special structure of our experimental setup. For any 
certificate, we are only interested in the roots between the 
current time and the next bending event, when the certificate 
becomes invalid. Many certificates do not have roots in this 
interval, but may have roots outside. The current implemen- 
tation first computes all real roots and thereafter discards the 
ones that lie outside the mentioned interval. 

Optimization 3. We use Descartes' rule of sign to cer- 
tify the non-existence of roots in the interval until the next 
bending event, and if successful, we skip the root isolation 
algorithm. 

The fifth column of Table [4] shows the improvement. More 
than ninety percent of the certificates that do not have a root 
before the next bending event are filtered out. As a final 
improvement, we avoid isolating the roots of the same poly- 
nomial multiple times. 

Optimization 4. We store polynomials together with 
their real roots in the interval until the next bending event 
in cache, which is cleared at the next bending event. 

We see in the sixth column of Table [4] that the cache yields 
another substantial speed-up, which suggests that certificates 
are frequently devalidated and revalidated during the run- 
time of the algorithm. We remark that also the kds kernel 
would benefit from caching. Comparing the running times 
for maintaining the alpha complex before and after the four 
steps of optimization, we see that the performance improves 
by roughly a factor of 15. Moreover, compared to maintain- 
ing the Delaunay triangulation, the optimized algorithm is 
slower by a factor up to 4. It is no surprise that the extension 



to alpha complexes is expensive. After all, it requires ad- 
ditional radius certificates, which have higher degrees than 
the flip certificates needed to maintain the Delaunay trian- 
gulation. We have demonstrated that with some algorithmic 
engineering, the overhead needed for alpha complexes can 
be kept within a moderate bound. 

6 Discussion 

The main contributions of this paper are a kinetic algorithm 
for alpha complexes, its use to construct a space-time repre- 
sentation of a spatial sorting process — called a medusa — 
and the implementation of a moderately fast but correct soft- 
ware. In a companion paper 0, we have demonstrated that 
the theory of persistent homology applied to the time func- 
tion on the medusa quantifies the sorting process by measur- 
ing its topological features. While the work in this and the 
companion papers is limited to simulated data, an applica- 
tion of our methods to observed biological data is under way. 
There is no theoretical obstacle to generalizing our algorithm 
and its implementation to the weighted case, in which dif- 
ferent Voronoi regions are restricted to within different size 
balls. A more challenging problem is the restriction to within 
bodies different than balls, e.g. arbitrarily oriented ellipsoids. 

The first author of this paper made a considerable effort to 
accelerate the implementation of the kinetic alpha complex 
algorithm, since this was necessary to compute examples of 
reasonable size in acceptable time; the instances computed in 
each took about 4 hours with our best configuration. Nev- 
ertheless, there are opportunities to further speed up the soft- 
ware, in particular on the level of the algebraic kernel. For 
example, it would be desirable to restrict the root isolation 
method to within a given interval, without wasting any time 
on roots outside this interval. Similarly, we expect an im- 
provement from implementing the filter for ruling out empty 
intervals in certified approximate arithmetic. We believe that 
kinetic data structures are an important tool in the topologi- 
cal analysis of time-varying shapes. We hope that our work 
on cell segregation initiates further work on such data. To 
facilitate this research, it would be useful if our extension 
of Coal's package on kinetic data structures is transformed 
from an experimental branch to a redesign of the package. It 
is desirable that such a redesign solves the problem of degen- 
eracies in the implementation of kinetic Delaunay triangula- 
tions and alpha complexes. Except for some special cases, 
the current versions of both algorithms are not guaranteed to 
work correctly when two or more events happen at exactly 
the same time. 

Acknowledgments 

The authors thank Viktoriia Sharmanska for discussions and help 
with a 2-dimensional prototype of our implementation. 



9 



References 



[1] E. BERBERICH, M. HEMMER AND M. KERBER. A generic alge- 
braic kernel for non-linear geometric applications. In "Proc. 27th 
Ann. Sympos. Comput. Geom., 2010", 179-186. 

[2] M. de Berg, 0. Cheong, M. van Krefeld and M. Over- 
mars. Computational Geometry: Algorithms and Applications. 
Springer- Verlag, Berlin, Germany, 2008. 

[3] H. EDELSBRUNNER. Geometry and Topology for Mesh Generation. 
Cambridge Univ. Press, Cambridge, England, 2001. 

[4] H. EDELSBRUNNER AND J. L. HARER. Computational Topology. 
An Introduction. Amer. Math. Soc, Providence, Rhode Island, 2010. 

[5] H. EDELSBRUNNER, C.-P. HEISENBERG, M. KERBER AND G. 
KRENS. The medusa of spatial sorting: topological construction. 
larXiv:1207.6474l 2012. 

[6] H. EDELSBRUNNER AND E. MUCKE. Three-dimensional alpha 
shapes. ACM Trans. Graphics 13 (1994), 43-72. 

[7] A. ElGENWILLIG, V. SHARMA AND C. YAP. Almost tight recursion 
tree bounds for the Descartes method. In "Internat. Sympos. Symbol. 
Alg. Comput., 2006", 71-78. 

[8] C.-P. HEISENBERG AND G. KRENS. Cell sorting in development. 
Current Topics in Developmental Biology, to appear. 

[9] M. HEMMER, E. TSIGARIDAS, Z. ZAFEIRAKOPOULOS, I. EMIRIS, 
M. KARAVELAS AND B. MOURRAIN. Experimental evaluation and 
cross-benchmarking of univariate real solvers. In "Proc. Conf. Sym- 
bol. Numeric Comput., 2009", 45-54. 

[10] L. Kettner, K. Mehlhorn, S. Pion, S. Schirra and C. Yap. 
Classroom examples of robustness problems in geometric computa- 
tions. Computational Geometry: Theory and Applications 40 (2008), 
61-78. 

[11] D. RUSSEL. Kinetic Data Structures in Practice. PhD Thesis, Stan- 
ford University, California, 2007. 

[12] D. RUSSEL. Kinetic Data Structures. CGAL User and Reference 
Manual, edition 4.0, 2012. 

[13] E. WEISSTEIN. Circumsphere. 

mathworld . wolfram . com/ Circumsphere . html ac- 
cessed, 2011. 

[14] WlKIPEDIA, THE FREE ENCYCLOPEDIA. Circumscribed sphere. 

[en . wikipedia . org/wiki/Circumscribed_circle#Circumcircle_equations| 
accessed September 11, 2012. 

[15] C. Yap. Towards exact geometric computation. Computational Ge- 
ometry: Theory and Applications 7 (1997), 3-23. 

[16] C. Yap. Fundamental Problems of Algorithmic Algebra. Oxford Uni- 
versity Press, Oxford, England, 2000. 



10 



